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ABSTRACT 

An Eulerian TVD code and a Lagrangian SPH code are used to simulate the 
off-axis collision of equal-mass main sequence stars in order to address the question 
of whether stellar mergers can produce a remnant star where the interior has been 
replenished with hydrogen due to significant mixing. Each parent main sequence star 
is chosen to be found near the turnoff, with hydrogen depleted in the core, and is 
modelled with a M — 0.8M Q realistic stellar model and as a n = 3 polytrope. An 
ideal fluid description with adiabatic index 7 = 5/3 is used for all hydrodynamic 
calculations. We found good agreement between the simulations for the polytropic 
case, with the remnant showing strong, non-local mixing throughout. In the interior 
quarter of the mass, ~ 35% is mixed in from larger radii and on average the remnant 
is ~ 50% fully mixed. For the realistic model, we found less mixing, particularly in 
the interior and in the SPH simulation. In the inner quarter, ~ 20% of the contained 
mass in the TVD case, but only ~ 3% in the SPH one is mixed in from outside. The 
simulations give consistent results for the overall profile of the merger remnant and the 
amount of mass loss, but the differences in mixing suggests that the intrinsic difference 
between grid and particle based schemes remains a possible artifact. We conclude that 
both the TVD and SPH schemes can be used equally well for problems that are best 
suited to their strengths and that care should be taken in interpreting results about 
fluid mixing. 

Key words: blue stragglers - globular clusters: general - hydrodynamics - methods: 
numerical - stars: evolution - stellar dynamics 



1 INTRODUCTION 

In dense stellar systems, such as globular clusters, galac- 
tic nuclei and star forming regions, direct collisions between 
stars occur quite frequently. These collisions can modify the 
stellar populations of the system by creating objects (e.g. 
blue stragglers, c ataclysmic variables and millisecond pul- 
sars; [Bailyn 1995), or by de stroying objects (e.g. bright gi- 
ants; iBailev fc Davieslll999h . These kinds of strong interac- 
tions between stars can also modify the overall evolution 
of the system by changing the dynamics of the system and 
its energy budget. The study of stellar collisio ns has be- 
come critical to the study of dense stellar systems l|Sills et al\ 
l2003h . 
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In order to understand how collisions can modify stellar 
populations and dynamics, we need to understand the colli- 
sions themselves. We need to understand when collisions will 
destroy the two stars involved, and when a new star will be 
created. If a star is created, how much mass has been lost 
from the system? What is the composition and structure of 
the new star? The answers to these questions depend on the 
two stars involved in the collision, their impact parameter, 
and the velocity of the collision. 

Consider for example, blue straggler stars, which are 
probable stellar collision products. These are stars in clusters 
that are bluer, brighter, and more massive than the main 
sequence turnoff, and hence must be rejuvenated stars. We 
see blue stragglers in all globular clusters that have been 
observed carefully, and they are also found in many open 
clusters and other environments. The number is typically 
a few 10's per cluster, but can get as high as a few hun- 
dred. Based on their numbers, these stars must have main 
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sequence lifetimes on the order of Gyrs, which implies a sig- 
nificant amount of hydrogen in their cores. Is that hydrogen 
mixed into the core during the stellar collision, or is it there 
already? Only a detailed understand of the hydrodynami- 
cal mixing processes can answer that question. Modelling 
the collision itself is a hydrodynamical problem because of 
the timescales involved. The merger timescale is typically on 
the order of hours, over which time the fluid is very close to 
ideal, with no radiative or diffusive effects of importance. 

The first hydrodynamic simulations of stellar collisions 
were done with low-reso lution grid codes, in one or two 
dimensions |Mathisl 1967). These earliest simulations were 
head-on collisions between equal-mass main sequence stars 
at a variety of impact velocities. Attention quickly turned to 
collisions between compact objects ( white dwarfs or neutron 
stars ) and main sequence stars (e.g. IShara fc Shavivlll977l. 



1978), and collisions involving giant branch stars (iTuchmanl 



1985). Grid methods are used today mostly for collisions in- 



volving neutron stars, where general relativity needs to be 
considered (e.g. ICalder fc Wandl2002l ). 

Smoothed part icle hydrody nami cs (SPH) was invented 
indep endently by iLucvl (1977) and iGingold fc Monaghanl 
l|l977h . The first SPH simulations o f stella r collisions were 
published in 1987 by iBenz fc Hills! (|l987l ). who looked at 
a variety of collisions involving equal-mass main sequence 
stars. Since then, the SPH method has been applied exten- 
sively to a wide range of collision scenarios. A very thor- 
ough revie w of collisions involvin g main sequence stars can 
be found in iFreitag fc Bend (120051 ). Almost all possible com- 
binations of stel lar types have been co llided, from pre-main 
sequence stars llLavcock fc Sills! 120051 ) through giants (e.g. 



iRasio fc Shapiro 



1991J) and even black holes and quark stars 



( Kluzniak fc Leell2002l ). 

The SPH method has some intrinsic advantages for 
studying stellar collisions. The main advantage is that the 
SPH particles follow the fluid flow. Therefore, all the compu- 
tational effort is spent on the regions of interest rather than 
in calculating empty grid cells. Some stellar collisions, par- 
ticularly those with large impact parameters and low rela- 
tive velocities, result in a very slow in-spiral of the two stars, 
with many periastron passages and large apastron distances. 
A grid code typically has to sacrifice resolution for coverage 
of the large physical space traversed by the stars involved in 
the collision. 

However, the SPH method also has some limitations 
which could affect our understanding of the outcome of the 
collision. Shocks are often not well resolved, in particular the 
weak ones which are expected in stellar mergers. Artificial 
viscosity is implemented to prevent unphysical oscillations, 
but it can cause spurious exchange of angular momentum 
between shear layers in a flow. A careful choice of the for- 
mulation of the artificial visocity is needed to address these 
concerns. In addition, there is usually no treatment of tur- 
bulent mixing. Since the final structure and subsequent evo- 
lution of a collision product can depend quite sensitively on 
its composition, getting the amount of fluid mixing right is 
very important. 

In grid codes, cells numerically force mixing of all mat- 
ter that enters. This will result in overmixing. On the other 
end of the scale, in SPH codes, particles have a definite iden- 
tity and entropy. Mixing in real fluids occurs through a tur- 
bulent cascade, where at the end some fluid elements gain 



entropy and some lose entropy. This process is not permit- 
ted in SPH, except in a coarse grained sense. In the presence 
of gravity, the residual identity could lead to sedimentation, 
which can undo the coarse grained mixing. From a physical 
perspective, we expect grid codes to overmix, and SPH to 
undermix, with the true answer in between. 

Some comparisons between SPH and grid cod es have 
been done for stellar collisions l|Davies et al.l [19931 ). They 
collided a low-mass unevolved main sequence star (modelled 
as a n=3/2 polytrope) with a white dwarf. They found that 
the global values of the simulations were in good agreement, 
but that some local values (e.g. the density distribution near 
the white dwarf) were slightly different. They also found 
that the largest differences between their simulations was 
caused by different treatments of the potential well of the 
white dwarf, which was treated as a point mass in both sim- 
ulations. The treatment of entropy and its effect on mixing 
was not investigated in that paper. 

Comparisons between SPH and grid codes for other 
hydrod ynamic problem s have also been published. For ex- 
ample, Bate fc Burkcrt ( 1997) looked at resolution require- 
ments for_j2otfi_meth£ds_jn_amol^^ cloud collapse calcu- 
lation iDi G. Sigalotti fc Klapd l|l997l ). They compare both 
first- and second-order hydrodynamic simulations of a col- 
lapse of a protostellar binary system with SPH, and find that 
the second-order codes produce results very similar to the 
SPH results. The main differences occur in the treatment 
of rotation. Accretio n flows were studied with both SPH 
and TVD methods bv lMolteni. Ryu, fc Chakrabartil (|l996h . 
who found that the general characteristics of the flow were 
the same, although the details were slightly different. They 
found that the location of the shocks were probably better 
obtained in SPH, but that the shocks themselves were better 
modelled using TVD. Comparisons between SPH and Eule- 
rian codes were a lso applied to dynamical bar i nstab ilities in 
a rotating star bv lSmith. Houser. fc Centrella ( 1996), and i n 
a variety of cosmological simulations bv lKang et al.l (| 1994h . 
The overall result of these studies is that both SPH and grid 
codes provide accurate descriptions of the hydrodynamics of 
most astrophysical situations. However, different scenarios 
are better treated with different methods, depending on the 
science that one wishes to do with the simulation. 

In this paper, we compare the results of stellar collision 
calculations perfo rmed using an Eulerian TVD grid code 
ilTrac fc Penll2003l ) and a Lagrangian SPH code (|Sills et al\ 
2002). We have simulated the off-axis collision of two equal- 
mass main sequence stars with Mo = 0.8Mq and Rq = 
0.955i?Q. Each parent main sequence star is chosen to be 
found near the turnoff, with hydrogen depleted in the core, 
and is modelled with a M = 0.8Mq realistic stellar model 
and as a n = 3 polytrope. We concentrate on a comparison 
of the mixing characteristics in the collision products in or- 
der to address the question of whether stellar mergers can 
produce a remnant star where the interior has been replen- 
ished with hydrogen due to significant mixing. In sections 
2 and 3 we describe the two hydrodynamic methods and in 
section 4 the numerical setup. Section 5 gives the compar- 
isons and we draw conclusions from those results in section 
6. 
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2 THE EULERIAN TVD CODE 

A three-dimensio nal Eulerian hydrodynamic code 
l|Trac fc Pen! | 20031 ) is used to compute the flux of the 
conserved fluid variables of mass, momentum, and total 
energy across a Cartesian grid in discrete time steps. For 
simplicity of notation, we can write the Euler equations in 
vector form as 

du dF(u) 



+ 



S(u) 



(1) 



dt dx 

where u contains the conserved fluid variables per unit vol- 
ume, F holds the hydrodynamic flux terms, and S repre- 
sents the gravitational source terms. The conserved cell- 
averaged quantities u„ = u(x n ) are defined at integer grid 
cell centres x n — k) while the fluxes F n+1 / 2 are defined 
at cell boundaries. The TVD hy dro code uses an operator 
splitting technique l|Strandll968T ) to separately solve the hy- 
drodynamic and gravitational terms. Furthermore, the hy- 
drodynamic step is dimensionally split such that flux updat- 
ing is done in one dimension at a time. 

In the hydrodynamic step, time integration is obtained 
using a second-order Runge-Kutta scheme, where a half-step 
is first performed, 
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and then followed by a full-step, 
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The second-order accurate fluxes F i+1 / 2 at cell boundaries 
are calcul ated using the total variation diminishing (TVD) 
condition (|Hartenlll983l ). The hydrodynamic step is conser- 
vative in that flux taken out of one cell is added to the neigh- 
bouring cells and the conservation of mass, momentum, and 
energy therefore ensures the correct shock propagation. 

To address the question of chemical mixing, the hy- 
dro code is implemented along with a particle-mesh scheme 
where tests particles of fixed mass flowing along the fluid 
velocity field lines can be used to track passively ad- 
vected quantities like chemical composition. The particle- 
mesh scheme is coupled to the Eulerian TVD scheme in the 
following manner. The hydrodynamic scheme is dimension- 
ally split and in order to synchronize the particles with the 
fluid, the particle-mesh scheme must also be dimensionally 
split. Consider solving for the fluid flow in the x direction. 
While the fluid velocity at cell centers can be defined by 
dividing the cell-averaged momentum by the cell mass, a 
more robust measure is that of the velocity at which mass 
is flowing across cell boundaries. The fluid velocity at cell 
boundaries can be calculated as 



v x, i+1/2 



i+1/2 



(4) 



where the fluxes F = F(p) come from solving the mass conti- 
nuity equation. Both the half-step velocity field v(xY x i+1 ; 2 

and full-step velocity field v(x) t J~^ t y 2 are nee ded for the 
particle-mesh scheme. The particles are then advected in 
the x direction using a second-order Runge-Kutta time in- 
tegration with a half step, 

i+At/2 t . , t 

X ' = X + V(X 



followed by a full step, 

x *+** = x * + v ( x t+At/2 ) x +At/2 At. (6) 

The velocity of the particle at position x is interpolated from 
the fluid velocity field using a 'cloud-in-cell' (CIC) technique 
(jHocknev &i Eastwoodllf988l ). The CIC interpolation is also 
used to determine local quantities such as density, pressure, 
entropy, and gravitational potential associated with each 
particle. The coupling of the particle-mesh algorithm to the 
Eulerian TVD code has the benefit of being able to track 
fluid variables like in an SPH code, while still being able to 
solve the fluid equations using the Eulerian approach. 



3 THE LAGRANGIAN SPH CODE 

The SPH co de is the pa r allel v ersion of the code described 
in detail in iBate et all (|l995h . It was parallelized using 
OpenMP. Our three dimensional code uses a tree to solve for 
the gravitational f orces and to find the nearest neighbours 
ijBenz et al. lll99d V We use an adiabatic equation of state, 
and the thermodynamic quantities are evolved by following 
the change in internal energy using the following equation: 



dui 
~dt 



3=1 



va ■ ViWi 



(7) 



(5) 



where u is the internal energy of a particle, m is its mass, 
P is pressure, p is density, n is the artificial viscosity, v 
is velocity, and W is the SPH smoothing kernel. Both the 
smoothing length and the number of nei ghbours can c hange 
in time and space. We used the Balsara (Balsara 1995) form 
of artificial viscosity with a = f3 — 5/6, following the advice 
of lLombardi et al.1 (l999). We ran a number of simulations 
with other forms of the artificial viscosity and with different 
values for the coefficients of the Balsara form. As expected, 
the differences between the implementations of the artificial 
viscosity were slight, except for the completely unphysical 
case of no artificial viscosity. 

Standard SPH treatments which do not properly ac- 
count for the change in smoothing length as a function of 
time or space can either evolve energy or entropy properly, 
but not the other variable. Since our SPH code evolves the 
energy equation, we are likely to have errors in the determi- 
nation of entropy and therefore may have unphysical fluid 
mixing in our simulations. However, earlier simulations of 
stellar collisions done with this code (Sills et al. 2002) have 
been compared to simulations by Lombardi et al. (1997) per- 
formed with an SPH code that evolves the entropy equation. 
The resulting hydrogen abundance profiles of the collision 
products were almost identical, and so we feel that the dif- 
ference between the entropy and energy treatment of the 
thermodynamics is minimal in these circumstances. 



4 HYDRODYNAMIC SIMULATIONS 

We have simulated the off-axis collision of two equal-mass 
main sequence stars with Mo = 0.8Mq and Ro = 0.955-Rq. 
The parent stars are chosen to be found near the turnoff, 
with hydrogen depleted in their cores, in order to address the 
question of whether stellar mergers can produce a remnant 
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Figure 1. Initial mass profiles of a Mo = 0.8Mq and Rq = 
O.955.R0 main sequence star constructed with (a) a realistic stel- 
lar model (solid lines) and (b) as an n = 3 polytrope (dashed 
lines). The TVD (blue) initial conditions exactly match the input 
models and the SPH (red) initial conditions are also in very good 
agreement. The SPH particles have been randomly subsampled 
by a factor of 5 for clarity. The entropic variable A = P/p 1 , where 
7 = 5/3 is the ideal gas adiabatic index, is used in place of the 
entropy. Units are in cgs. 



star where the hydrogen core has been replenished from sig- 
nificant chemical mixing. The main sequence stars are con- 
structed using two different methods and the initial mass 
profiles are compared in Figure [T] First, a realistic model 
is calculated with th e Yale stellar evolution code (YREC; 
iGuenther et aZ.|[l992l ). The YREC stars are found near the 
turnoff with an age of 13.5 Gyr. In the central ~ 25% of the 
mass, formerly the hydrogen core, the hydrogen has been 



mostly converted to helium. Second, a polytropic model 
with index n = 3 gives a good approximation to main se- 
quence stars near the turnoff that have relatively little mass 
in their convective envelopes. Note that in contrast, n=3/2 
polytropes are good approximations of convective regions, 
and very good for fully convective or mostly convective un- 
evolved low mass stars. In both of our chosen stellar models, 
approximately 90% of the total mass is contained within 
r < 0.5-Ro- However, the realistic model is more centrally 
condensed with steeper gradients than the polytropic one. 

The parent stars are set up on zero-energy parabolic or- 
bits with a pericentre separation equal to 0.25i?o and the ini- 
tial trajectories are calculated assuming point masses. The 
orbital plane is taken to coincide with the x — y plane and 
the stars are initially separated by a distance of 3.75-Ro. The 
dynamical time for the stars, defined to be 

T dyn =—^=, (8) 



fGp 

where p is the average density, is approximately 1 physical 
hour and thus, modelling of the collision is a hydrodynamic 
problem. 

In the Eulerian simulations, the collision is computed on 
a high resolution Cartesian grid with 1024 3 cells. Initially, 
each parent star has a radius of Rq = 192 grid cells and is 
assigned N p = 256 3 test particles. The hydro grid has non- 
periodic boundaries which allow only one-way transmission 
and gas exiting the grid can not re-enter. It is important 
that the parent stars are chosen to be small enough that the 
mass loss is still accurately captured on the finite grid. The 
FFT Poisson solver can also be made non-periodic by using 
a gravity grid with dimensions that are twice as large as the 
hydro grid and treating the padded cells as vacuum with 
zero density. However, since the remnant star occupies only 
a fraction of the hydro grid, a more optimal approach is to 
make the Poisson solver only pseudo non-periodic and use 
a gravity grid with dimensions that are 1.5 times as large 
as the hydro grid. We impose a minimum density of 10~ 6 
g/cm 3 whereby gas cells with lower density are considered 
part of the vacuum and all fluxes are set to zero. This density 
threshold is found to reduce spuri ous shock heating as the 
stars are advected across the grid (|Trac fc Penll2003r ) and is 
still a few orders of magnitude smaller than densities of in- 
terest. An ideal fluid approximation with an adiabatic index 
of 7 = 5/3 is used in the simulations. The simulations were 
run at NCSA0 on a shared-memory SGI Altix with Itanium 
2 processors. Each simulation used 64 processors and 64 GB 
of memory and required approximately 8000 cpu hours to 
complete. 

In the SPH simulations, we used approximately 50 
000 particles per star, or 100 000 in total. The particles 
were initially distributed on an equally-spaced grid, and 
their masses were varied until the density profile matched 
that of the stellar model we wished to use. The smoothing 
length is varied to keep the number of neighbours approx- 
imately constant at about 50. The particle masses range 
from 3.2 x 10~ 2 Af Q to 7.8 x 1O" 8 M with the highest mass 
particles at the centre of the original star. The simulations 
were run using SHARCNET0 facilities on a 4-CPU shared- 



1 http://www.ncsa.uiuc.edu 

2 http://www.sharcnet.ca 
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Figure 2. Density plots of the merger remnants for the realistic 
stellar model from the TVD (top) and SPH (bottom) simulations. 
Thin slices approximately 0.05-Rq thick are taken through the 
x — y and x — z mid-planes. The orbital plane coincides with the 
x — y mid-plane. Both simulations produce a single rotating oblate 
with a mass of M = 1.5Mq. 

memory node of a Compaq SC40 Alpha with 667 MHz pro- 
cessors. Each simulation used less than 2 GB of memory and 
required approximately 150 CPU hours. Note the substan- 
tially smaller resource requirements of the SPH method for 
this kind of simulation compared to the grid method. 



5 RESULTS 

Both TVD and SPH simulations of the off-axis collision be- 
tween Mo = 0.8M Q main sequence stars and between n = 3 
polytropes produced a single merger remnant with mass 
M* = 1.5M0that establishes hydrostatic equilibrium in ap- 
proximately lOrdyn- Note that from here on, M, refers to the 
mass of bound gas and does not include ejected gas. During 
the merging process, the outer envelopes of the parent stars 
are shock-heated and O.IM© of gas gets ejected into the in- 
terstellar medium (ISM) . The resulting remnant is a rotating 
oblate that shows no indication of having a disc. The TVD 
and SPH simulations find similar mass loss with each other 
and also for the two stellar models. The latter is due to the 
fact that both stellar models have similar structural profiles 
at larger radii. We do not attempt to quantify it to higher 
precision because of the difficulty in precisely determining 
what gas is bound to the remnant. We show some images 
of the remnant in Figure [5] taken at time t — 10r c i yn . All 
future figures show results also taken at this final time. 

5.1 Mass profiles 

Figures [3^, and[3p show the mass profiles of the merger rem- 
nant from the collision between main sequence stars and 
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Figure 3. Mass profiles of the merger remnant from the off-axis 
collision between equal mass Mo = 0.8Mq main sequence stars 
(a) and between n = 3 polytropes (b), simulated using a TVD 
(blue) and a SPH (red) code. Note that ideal gas adiabatic index 
of 7 = 5/3 is used to define the entropic variable A = Pip 1 . Units 
are in cgs. 

between polytropes, respectively. The two lines span the 1- 
sigma spread in the quantities in each mass bin. In this and 
all future figures, the red lines give the results of the SPH 
simulations, while the blue lines are the TVD results. 

The profiles are constructed by first rank ordering the 
grid cells or particles in terms of the effective potential and 
then binning the structural variables. The entropic variable, 

A = P/p\ (9) 

where 7 = 5/3 is the ideal gas adiabatic index, is used 
in place of the entropy. We also check if the particle-mesh 
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scheme in the Eulerian code is working accurately by mea- 
suring the mass profiles for the test particles. The CIC in- 
terpolation technique is used to assign local gas quantities 
to the particles and if the particles are indeed tracing the 
gas, then the mass profiles should be the same. The particle- 
mesh algorithm does an accurate job of advecting the test 
particles along the velocity field lines as we find agreement 
at the percent level over the entire mass range. 

The TVD and SPH results for the polytropic case are 
in good agreement at all mass scales. The remnant star is 
a rotating oblate with a mass of M* = 1.5M Q and is it- 
self not a polytrope. The central density and pressure have 
decreased and the entropic variable floor has risen, relative 
to the corresponding values in the parent polytropes. These 
changes can come from a combination of relaxation, shock 
heating, and mixing. The initial polytropes and the final 
remnant all have structural profiles that are smooth and 
featureless, and the modest demand on dynamic range is 
more than accounted for by the high resolution simulations. 
Numerical relaxation in the core has been f ound to be o nly 
a few percent at these high resolutions. In I Trad l|2004 ). a 
n — 3 polytrope with a radius of R — 96 grid cells, which is 
lower resolution than that used in this paper, was advected 
diagonally across the x — y orbital plane in over 1600 times 
steps to test the TVD code. For the inner 5% of the en- 
closed mass, the decrease in density and pressure and the 
increase in entropy was found to be < 5%, while for the rest 
of the polytrope, negligible change was found. Similarly, a 
test of the SPH code, where an isolated n — 3 polytrope 
was evolved for 10 dynamical times, showed that the mix- 
ing in the inner 25% of the enclosed mass was < 2%. Since 
shock heating is expected to be minimal in the core because 
of the low Mach numbers, the change in entropy suggests 
that some mixing may have taken place. We will address 
the question of mixing in more details in section 15.21 

For the main sequence stars, both codes give similar re- 
sults for the outer two-thirds of the total mass, but there 
are interesting differences in the interior. Like in the poly- 
tropic case, the remnant star is a rotating oblate with a mass 
of M* = 1.5Af© and some structure is retained from that of 
the parent stars. In the TVD results, the central density and 
pressure are lower and the entropic variable floor is higher 
compared to the corresponding SPH values. These observed 
differences are consistent with the differences between the 
Eulerian and Lagrangian approaches of the two simulations. 
In the Eulerian case, there will always be some numerical dif- 
fusion between adjacent grid cells and while this is generally 
alleviated by increasing the resolution, it is difficult to fully 
remove this effect here. In the hydrogen-depleted core of the 
main sequence stars, the structural variables change rapidly 
with mass or radius and the demand on dynamic range is 
challenging for a grid code without refinement meshes. In 
the SPH results, it appears that the cores of the parent stars 
have settled to the bottom of the potential well of the rem- 
nant star with little change in entropy. This is consistent, 
but not necessarily confirmation of the argument that since 
the weak shocks present in these type of mergers are not 
well-resolved in SPH simulations, particles will simply follow 
lines of constant entropy and may experience sedimentation. 

The particular implementation of the SPH equations 
used in this paper was based on non-equal-mass particles. 
For this reason, the inner regions of the parent stars consist 



of a few massive particles, and the density distribution re- 
quires that the particles in the inner regions of the realistic 
stars are even more massive than those in the polytrope. As 
a result, when the particles are binned in the collision prod- 
uct, the innermost bins have very few particles, resulting in 
the 'teeth' seen in the top half of figure [3] The resolution in 
the inner regions is somewhat worse than that of the grid 
code simulations. We should point out, however, that the 
overall results of the SPH simulations are not significantly 
affected by this lack of resolution. SPH simulations that use 
equal-mass particles show the same behaviour as the ones 
presented here, in that the cores of the parent stars do not 
change their entropy and si nk to the centre of the product 
|Lombardi et al.lll996l . l2002h . 

5.2 Mass mixing 

One of the advantages of Lagrangian SPH simulations is 
that the particles of constant mass are directly tracked al- 
lowing one to probe the mixing of mass and the change in 
entropy, and to determine the final chemical composition of 
the merger remnant. This is now possible for Eulerian hydro 
codes coupled with a particle-mesh algorithm. 

The mixing of mass is closely examined with the fol- 
lowing procedure. The particles are ranked ordered using 
their effective potential and each particle is assigned a mass 
coordinate based on its ranking and this mass coordinate is 
simply the value of the enclosed mass. In Figure|4]the initial 
and final mass coordinates are compared for both the TVD 
test particles and the SPH particles. The particles are binned 
in terms of the final mass M/ in order to see where the mass 
in a given region of the remnant is coming from. For each 
bin the 50% value of M; is found and the median curve is 
illustrated by the solid line. Similarly, the long-dashed lines 
enclose 50% of the scatter, while the short-dashed lines en- 
close 90% of the scatter. 

Both simulations of the merger of n — 3 polytropes 
produced substantial mixing. For the inner region of the 
remnant given by Mf/M» < 0.25 where the replenishing of 
hydrogen is important, the median curve for Mi deviates 
from the one-to-one line as mass is mixed into the interior, 
mainly from larger radii as the mixing from below has a 
hard limit. Considering this region as a whole, 38% of the 
contained mass in the TVD case and 32% in the SPH one 
is mixed in from outside. The average amount of mixing for 
the entire remnant can be estimated using the average width 
AM, of the zone enclosing half of the scatter, as bounded 
by the long-dashed lines. For complete mixing, the width is 
given by AM, = M»/2. In the TVD case, the remnant is 
estimated to be ~ 56% fully mixed while for the SPH case, 
we find that it is ~ 46% fully mixed. 

The merger of Mo = 0.8Mq main sequence stars results 
in less mixing, particularly in the interior of the remnant 
and in the SPH simulation. Considering the interior region 
Mf/M, < 0.25 as a whole, 21% of the contained mass in 
the TVD case, but only 3% in the SPH one is mixed in from 
outside. On average, the TVD remnant is estimated to be 
~ 48% fully mixed while the SPH remnant is only ~ 27% 
fully mixed. The weaker mixing seen here in this case com- 
pared to the polytrope merger is due to the structural dif- 
ferences between the two stellar models. The realistic main 
sequence star is more centrally condensed and it is relatively 
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Figure 4. Mass mixing found in the mergers of main sequence 
stars (a) and polytropes (b). The median (solid), 50% scatter 
(long-dashed), and 90% scatter (short-dashed) values are plotted. 

harder to penetrate the hydrogen-depleted core. During the 
merging process, the dense interior of each parent star acts 
like a semi-compact object and remains largely intact until 
the two cores reach the bottom of the potential well and 
merge with one another. The differences in mixing between 
the reali stic and polytropic mer gers have previously been 
found bv lSills fc Lombardil (ll997T ) for SPH simulations. The 
stronger mixing seen in the Eulerian case is consistent with 
the tendency of grid codes to overmix, whereas particle codes 
tend to undermix. The effect of this mixing on the hydrogen 
prof ile is discussed later in section 15.41 

lLombardi et al.l (Il996l ) have also studied the mixing 
of mass in SPH simulations of stellar collisions for vari- 
ous stellar models. For comparison, their closest model to 



ours involves the collision between n = 3 polytropes with 
Mo = IMq and Ro ~ IRq. When considering the interior 
quarter of the enclosed initial mass Mi /Mo, they found a 
mixing value of ~ 12%. Unlike our simulations, they used 
equal-mass particles, and so had approximately 7000 parti- 
cles in the inner 25% of the star. In our simulation, we have 
just over 1500 particles in that region of the merger rem- 
na nt, despite having more particles in our simulations over- 
all. HHE 

et al. (2002) show that lower resolution simulations 
result in more mixing (see their figure 1), simply because 
the regions of interest are by necessity more smoothed. We 
conclude that our results are in good general agreement with 
other SPH simulations of stellar collisions. 

5.3 Entropy change 

In Figure [5] the initial and final values of the entropic vari- 
able A = P/p"< of the TVD test particles and the SPH par- 
ticles are compared. Note that ideal gas adiabatic index of 
7 = 5/3 is used to define the entropic variable. In the ab- 
sence of shocks, mixing, and relaxation, the fluid elements 
in the parent stars should be simply sorted by entropy in 
the remnant star, with Af = Ai. 

In the polytropic merger, it is not surprising that signif- 
icant change in the entropic variable is found considering the 
amount of mixing that has occurred. The two simulations are 
also in good agreement, though the small differences suggest 
that the TVD remnant is more strongly mixed, as seen in 
the previous section, and more strongly shock heated, which 
is expected because grid codes are known to provide better 
capturing of shocks. Note that while strong mixing has been 
found to occur in the interior of the remnant, it has only 
raised the entropy floor by a small amount. 

With main sequence stars, the TVD results show an 
artificial increase in the entropy floor that is caused by nu- 
merical relaxation rather than mixing or shock heating. In 
comparison, the lowest entropy SPH particles show only rel- 
atively small changes with small dispersion. Overall, there is 
less change in the entropic variable in the realistic merger, 
but it is difficult to tell from the plots how much of the en- 
tropy increase is due to mixing and how much to shock heat- 
ing. To account for the differences between Figure [S^- and 
Figure [5p, mixing must play a substantial role in changing 
the entropy profile or there must be significantly more shock 
heating in the polytropic merger. Note that for some TVD 
test particles where shocks have not significantly raised the 
entropic variable, the decrease in local entropy due to mixing 
is observable. 

5.4 Hydrogen profile 

The hydrogen mass fraction X in the merger remnant from 
the collision between main sequence stars is plotted as a 
function of the enclosed mass in Figure [U] For compar- 
ison, the short-dashed line represents the hydrogen pro- 
file expected in the absence of mixing. The SPH remnant 
has a hydrogen profile that closely resembles the no mix- 
ing case. The hydrogen profile in the hydrogen-depleted 
interior is sensitive to mixing, but only weak mixing and 
small changes of < 5% to the hydrogen fraction are found. 
While stronger mixing is found at larger radii, it occurs be- 
tween regions with approximately similar initial hydrogen 
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Figure 5. Change in the entropic variable found in the mergers of 
main sequence stars (a) and polytropes (b). The median (solid), 
50% scatter (long-dashed), and 90% scatter (short-dashed) val- 
ues are plotted. The diagonal dotted line represents no change 
while the horizontal dotted line gives the entropic variable for the 
interior quarter of the remnant mass Mf/M* = 0.25. Note that 
ideal gas adiabatic index of 7 = 5/3 is used to define the entropic 
variable A = P/p 1 ■ Units are in cgs. 



abundance and this results in only small changes. In the 
TVD remnant the stronger, non-local mixing has noticeably 
changed the hydrogen profile. Considering the interior re- 
gion Mf /M* < 0.25 as a whole, the hydrogen fraction has 
on average increased by 24% and this is consistent with the 
previous finding where 21% of the mass is mixed in from the 
outside. 

We also include the results for the polytropic collision 
to highlight the difference in mixing, but since this is not 



Figure 6. Hydrogen mass fraction in the merger remnant from 
the collision between main sequence stars (solid) and between 
polytropes (long-dashed). The short-dashed black line represents 
the hydrogen profile expected in the absence of mixing. 

a realistic model, the initial chemical composition of the 
polytropes was taken from that of the Mo = 0.8Mq main 
sequence star. In this case, both simulations show stronger, 
non-local mixing and this has resulted in more changes to 
the hydrogen mass fraction throughout. 



6 DISCUSSION AND CONCLUSIONS 

The biggest difference we see between the SPH and TVD 
simulations is in the amount of mixing in the interior of the 
remnant star. In order for a collision product to be observed 
as a blue straggler, it must remain as a main sequence star 
for a reasonably long time. The SPH simulations predict 
a very small amount of mixing, and we know f rom stellar 
evolu tionary models of these collision products (Sill s et all 
119971 ) that the lifetimes of these collision products are quite 
short (~ 10 6 years). The TVD simulation of the same colli- 
sion produces a collision product which has a central hydro- 
gen abundance of X ~ 0.3, which will have a significantly 
longer lifetime (~ 10 s years). It is not clear from this work 
which answer is more correct. We probably need to turn to 
the observations of blue stragglers in a variety of clusters. 
We know already that simple blue straggler population syn- 
thesis models (based on results from SPH simulations) have 
some trouble reproducing the bri ght end of the observed dis- 
tribution in some clusters (ISills fc Bailvnlll999l) . but we also 
know that Eulerian schemes contain some numerical diffu- 
sion and hence over-predict the mixing. The true result is 
probably somewhere between these two cases. 

Note that the mixing effect will be strongest in the col- 
lision that we have chosen to model: a collision between two 
stars at the main sequence turnoff. This is, of course, why 
we chose to compare this particular collision, but we should 
caution the reader not to apply these results uniformly to 
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all collisions. Stars that are less evolved and that have flat- 
ter composition profiles will show the effect of any mixing 
less when they collide, since any material that is 'mixed' will 
have the same composition. The particular choice of hydro- 
dynamical method should not produce such a different result 
as it does for the highly evolved case shown here. 

The differences in mixing seen in the merger of main 
sequence stars and the merger of polytropes are due to the 
structural differences between the two stellar models. The 
realistic main sequence star is more centrally condensed and 
it is relatively harder to penetrate the hydrogen-depleted 
interior. Therefore, one expects that the larger the density 
contrast in the parent star, the weaker the mixing for a 
fixed collision velocity and geometry. Given the differences in 
results, it is important to use realistic stellar models rather 
than appr oximations like polytrop es when simulating stellar 



tnan approximations like polytropc 
collisions (|Sills fc Lombardilll997h . 



The Eulerian simulations have the advantage of not re- 
quiring any artificial viscosity or the need to adjust any 
smoothing lengths and therefore, the numerical solutions to 
the fluid equations, in particular the capturing of shocks, 
are generally robust. However, they have the disadvantage 
of limited dynamic range due to the fact that the fluid equa- 
tions were solved on a single Cartesian grid. One conse- 
quence of this is that only a small collisional parameter space 
can be considered and another is that the simulations often 
suffer some level of numerical diffusion in regions with large 
density gradients. In principle, the problems can be rectified 
by running larger volume and higher resolution simulations, 
but in practice this is not feasible because of limited com- 
putational resources. One remedy is to use adap tive mesh 
refine ment ( AMR) schemes in codes like A RT dKravtsov 
1993), Enzo (|Q'Shea et aZ.l2005l ). and FLASH (|Frvxell et al. 
200(3 ). One caution we make is that it is important to un- 
derstand any systematic problems in the transport of fluid 
across refinement mesh boundaries when studying mixing 
during stellar mergers. 

The SPH method is the hydrodynamical method of 
choice for modelling stellar collisions, because of its adapt- 
ability to a wide range of impact parameters, velocities, 
stellar masses, etc. One does not need to know where the 
material will go before starting the simulation, since the 
SPH particles are free to travel wherever they would like 
in space. The method has its drawbacks in terms of treat- 
ment of shocks and spurious diffusion of particles caused 
by the numerical methods. In addition, turbulent mixing 
is generally inhibi ted by the addition of artificial viscosity 
(jDolag et aLl2005l ) . We have shown in this paper that for the 
most part, the two methods produce the same results. The 
structure of the resulting collision product is very similar: 
the density and pressure profiles of the models are almost 
identical. However, the amount of fluid mixing is somewhat 
different between the two methods, which can result in dif- 
fering stellar lifetimes for the collision product. 

We conclude that both Eulerian and Lagrangian imple- 
mentations of fluid dynamics can be used to study stellar 
collisions, and that each code should be applied in a way to 
maximize its strengths. Care should be taken when inter- 
preting results related to fluid mixing, but all other prop- 
erties of the system are determined correctly by both ap- 
proaches. 
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